﻿/* Copyright 2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Marcus Cuda
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 */
using System;
using dnAnalytics.Math;
using dnAnalytics.Properties;

namespace dnAnalytics.LinearAlgebra.Decomposition
{
    /// <summary>
    /// Creates an orthogonal matrix using the modified Gram-Schmidt method.
    /// </summary>
    /// <remarks>The actual decompostiion is not done until one of the class'
    /// methods is invoked.</remarks>
    public class GramSchmidt : IQR
    {
        private readonly int mColumns;
        private readonly int mRows;
        private bool mComputed;
        private double mDeterminant = 1;
        private bool mIsFullRank = true;
        private Matrix mMatrix;
        private Matrix mQ;
        private Matrix mR;


        /// <summary>
        /// Constructs the GramSchmidt object for the given matrix.
        /// </summary>
        /// <param name="matrix">The matrix to decompose.</param>
        /// <remarks>The actual decompostiion is not done until one of the class'
        /// methods is invoked.</remarks>
        public GramSchmidt(Matrix matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }
            if (matrix.Rows < matrix.Columns)
            {
                throw new InvalidMatrixOperationException(Resources.RowsLessThanColumns);
            }
            mMatrix = matrix.Clone();
            mRows = mMatrix.Rows;
            mColumns = mMatrix.Columns;
        }

        #region IQR Members

        /// <summary>
        /// Returns the orthogonal Q matrix.
        /// </summary>
        /// <returns>The orthogonal Q matrix.</returns>
        public Matrix Q()
        {
            Compute();
            return mQ.Clone();
        }

        /// <summary>
        /// Returns the upper triangular factor R.
        /// </summary>
        /// <returns>The upper triangular factor R.</returns>
        public Matrix R()
        {
            Compute();
            return mR.Clone();
        }

        /// <summary>
        /// Determines whether the matrix is full rank or not.
        /// </summary>
        /// <returns><b>true</b> if the matrix is full rank; otherwise <b>false</b>.</returns>         
        public bool IsFullRank()
        {
            Compute();
            return mIsFullRank;
        }

        ///<summary>Calculates the determinant (absolute value) of the matrix.</summary>
        ///<returns>The determinant of the matrix.</returns>
        ///<exception cref="MatrixNotSquareException">If the matrix is not square.</exception>
        public double Determinant()
        {
            Compute();
            return mDeterminant;
        }

        /// <summary>
        /// Solves a system of linear equations, <b>AX = B</b>.
        /// </summary>
        /// <param name="input">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
        /// <returns>
        /// The left hand side <see cref="Matrix"/>, <b>X</b>.
        /// </returns>
        public Matrix Solve(Matrix input)
        {
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }
            if (mRows != input.Rows)
            {
                throw new NotConformableException("input", Resources.ParameterNotConformable);
            }

            Compute();
            Matrix result = mR.CreateMatrix(mQ.Columns, input.Columns);
            Solve(input, result);
            return result;
        }

        /// <summary>
        /// Solves a system of linear equations, <b>AX = B</b>.
        /// </summary>
        /// <param name="input">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
        /// <param name="result">The left hand side <see cref="Matrix"/>, <b>X</b>.</param>
        public void Solve(Matrix input, Matrix result)
        {
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (input.Columns != result.Columns)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (mRows != input.Rows)
            {
                throw new NotConformableException("input", Resources.ParameterNotConformable);
            }

            if (mColumns != result.Rows)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }
            Compute();

            Matrix b = input.Clone();
            int bn = b.Columns;

            // Compute Y = transpose(Q)*B
            double[] column = new double[mRows];
            for (int j = 0; j < bn; j++)
            {
                for (int k = 0; k < mRows; k++)
                {
                    column[k] = b.ValueAt(k, j);
                }
                for (int i = 0; i < mColumns; i++)
                {
                    double s = 0;
                    for (int k = 0; k < mRows; k++)
                    {
                        s += mQ.ValueAt(k, i)*column[k];
                    }
                    b.ValueAt(i, j, s);
                }
            }

            // Solve R*X = Y;
            for (int k = mColumns - 1; k >= 0; k--)
            {
                for (int j = 0; j < bn; j++)
                {
                    b.ValueAt(k, j, b.ValueAt(k, j)/mR.ValueAt(k, k));
                }
                for (int i = 0; i < k; i++)
                {
                    for (int j = 0; j < bn; j++)
                    {
                        b.ValueAt(i, j, b.ValueAt(i, j) - b.ValueAt(k, j)*mR.ValueAt(i, k));
                    }
                }
            }

            for (int i = 0; i < mR.Columns; i++)
            {
                for (int j = 0; j < input.Columns; j++)
                {
                    result.ValueAt(i, j, b.ValueAt(i, j));
                }
            }
        }

        /// <summary>
        /// Solves a system of linear equations, <b>Ax = b</b>.
        /// </summary>
        /// <param name="input">The right hand side vector, <b>b</b>.</param>
        /// <returns>
        /// The left hand side <see cref="Vector"/>, <b>x</b>.
        /// </returns>
        public Vector Solve(Vector input)
        {
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }

            if (mRows != input.Count)
            {
                throw new NotConformableException("input", Resources.ParameterNotConformable);
            }

            Compute();

            Vector result = mQ.CreateVector(mR.Columns);
            Solve(input, result);
            return result;
        }

        /// <summary>
        /// Solves a system of linear equations, <b>Ax = b</b>.
        /// </summary>
        /// <param name="input">The right hand side vector, <b>b</b>.</param>
        /// <param name="result">The left hand side <see cref="Matrix"/>, <b>x</b>.</param>
        public void Solve(Vector input, Vector result)
        {
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (mRows != input.Count)
            {
                throw new NotConformableException("input", Resources.ParameterNotConformable);
            }

            if (mColumns != result.Count)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            Compute();

            Vector b = input.Clone();

            // Compute Y = transpose(Q)*B
            double[] column = new double[mRows];
            for (int k = 0; k < mRows; k++)
            {
                column[k] = b[k];
            }
            for (int i = 0; i < mColumns; i++)
            {
                double s = 0;
                for (int k = 0; k < mRows; k++)
                {
                    s += mQ.ValueAt(k, i)*column[k];
                }
                b[i] = s;
            }

            // Solve R*X = Y;
            for (int k = mColumns - 1; k >= 0; k--)
            {
                b[k] /= mR.ValueAt(k, k);
                for (int i = 0; i < k; i++)
                {
                    b[i] -= b[k]*mR.ValueAt(i, k);
                }
            }

            for (int i = 0; i < mR.Columns; i++)
            {
                result[i] = b[i];
            }
        }

        #endregion

        private void Compute()
        {
            if (!mComputed)
            {
                mQ = mMatrix.CreateMatrix(mRows, mColumns);
                mR = mMatrix.CreateMatrix(mColumns, mColumns);
                for (int k = 0; k < mColumns; k++)
                {
                    double norm = mMatrix.GetColumn(k).PNorm(2);
                    mR.ValueAt(k, k, norm);
                    for (int i = 0; i < mRows; i++)
                    {
                        mQ.ValueAt(i, k, mMatrix.ValueAt(i, k)/norm);
                    }
                    for (int j = k + 1; j < mColumns; j++)
                    {
                        double dot = mQ.GetColumn(k).DotProduct(mMatrix.GetColumn(j));
                        mR.ValueAt(k, j, dot);
                        for (int i = 0; i < mRows; i++)
                        {
                            double value = mMatrix.ValueAt(i, j) - mQ.ValueAt(i, k)*dot;
                            mMatrix.ValueAt(i, j, value);
                        }
                    }
                }
                mMatrix = null;

                for (int i = 0; i < mColumns; i++)
                {
                    mDeterminant *= mR.ValueAt(i, i);
                    if (Precision.EqualsWithinDecimalPlaces(System.Math.Abs(mR.ValueAt(i, i)), 0.0, Constants.DecimalsForComparison))
                    {
                        mIsFullRank = false;
                        mDeterminant = 0;
                        break;
                    }
                }
                mDeterminant = System.Math.Abs(mDeterminant);
                mComputed = true;
            }
        }
    }
}
